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Executive  summary 

The  scope  of  the  work  as  proposed  to  the  funding  agency  was  to  carry  out  numerical  calculations  of  high-speed 
mixing  flows  using  GASP,  a  computational  fluid  dynamics  (CFD)  code.  The  ultimate  goal  is  to  use  GASP  for  pre¬ 
diction  and  optimization  of  chemical  laser  flows.  To  do  it,  one  must  have  confidence  in  the  performance  of  the 
hydrodynamic  and  mixing  part  of  the  code. 

For  complex  numerical  simulations,  it  is  essential  to  perform  code  validation  -  a  quantitative  comparison  of  nu¬ 
merical  results  with  experimental  results  confirming  that  the  code  faithfully  reproduces  the  physics  of  the  real  problem. 
Two  detailed  validation  exercises  were  performed  with  GASP: 

•  Simulation  of  a  shock-accelerated  mixing  flow.  This  validation  exercise  compared  the  numerical  results 
produced  by  GASP  with  highly-resolved  experimental  data  on  flows  subject  to  Richtmyer-Meshkov  instability 
(RMI),  the  preferred  test  problem  for  quantitative  assessment  of  the  performance  of  CFD  codes  in  prediction  of 
the  properties  of  high-speed  mixing  flows. 

The  validation  problem  revealed  that  GASP  can  faithfully  reproduce  all  the  large-scale  quantitative  properties 
of  the  flow.  We  also  gained  important  additional  insights  into  the  strengths  and  limitations  of  GASP  (and  other 
numerical  codes)  in  prediction  of  disordered  small  scales  in  flows  transitioning  to  turbulence. 

•  General  jet-in-crossflow  problem.  The  problem  of  a  supersonic  jet  discharging  into  a  supersonic  crossflow 
presents  a  significant  computational  challenge,  yet  provides  a  valuable  validation  exercise  when  the  results  are 
compared  with  experimental  data. 

The  conditions  for  this  exercise  were  chosen  to  achieve  two  goals: 

-  Provide  comparison  with  best  experimental  data  available  in  open  literature. 

-  Test  the  code  in  a  generic  configuration  that  can  be  easily  modified  to  represent  a  specific  chemical  laser 
nozzle  injection  scheme. 

The  agreement  with  the  experimental  data  was  quite  good,  yet  subtle  features  of  the  flow  that  were  not  apparent 
in  the  experiment  were  revealed  by  the  simulation. 

The  results  of  our  work  are  as  follows. 

•  We  found  the  hydrodynamic  model  used  in  GASP  to  be  suitable  for  prediction  of  mixing  in  chemical  laser  flows. 

•  We  found  the  error  in  prediction  of  the  geometry  of  large-scale  features  (e.g.,  counter-rotating  vortex  pairs)  not 
to  exceed  10-15%  when  compared  with  experiment. 
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•  We  found  that  a  fully  three-dimensional  model  is  advantageous  when  the  goal  is  to  correctly  assess  mixing. 

•  We  determined  the  grid  parameters  and  the  appropriate  conditions  to  run  the  simulation  to  produce  credible 
results  in  the  range  of  flow  regimes  consistent  with  the  operation  of  a  chemical  oxygen-iodine  laser  (COIL). 
GASP  has  been  validated  by  comparison  with  experinental  data. 

The  detailed  report  that  follows  begins  with  an  overview  of  the  two  validation  problems  we  considered,  followed 
by  the  description  of  the  computational  set-up  and  the  discussion  of  the  results.  The  work  summarized  here  has  been 
presented  at  AIAA  [1]  and  APS  [2,  3]  conferences,  as  well  as  in  a  journal  paper  currently  under  review  [4]. 


1  Overview  of  validation  problems 

The  two  validation  problems  were  selected  to  provide  good  quantitative  comparison  with  experiment.  There  are 
relatively  few  available  data  sets  characterizing  high-speed  mixing  flow  with  any  level  of  quantitative  detail.  Our 
first  validation  problem,  shock  acceleration  of  a  heavy  gas  cylinder  immersed  in  lighter  gas,  was  chosen  because  of 
relatively  simple  initial  conditions,  which  are  nominally  two  dimensional,  making  it  possible  to  run  a  two-dimensional 
simulation  and  compare  the  flow  features.  The  second  problem  is  much  more  complex  geometrically,  and  it  features 
flow  conditions  that,  while  being  sufficiently  generic,  are  very  close  to  those  in  a  chemical  laser  -  supersonic  injection 
of  gas  into  a  supersonic  flow. 

In  a  chemical  laser,  hydrodynamic  mixing  of  fuel  and  oxidizer  is  essential  to  obtain  a  population  inversion  of  the 
excited  state  at  an  optimum  rate.  If  mixing  is  poor,  then  the  chemical  reactions  occur  too  slowly,  producing  few  excited 
atoms  in  the  lasing  cavity.  On  the  other  hand,  if  the  mixing  is  too  fast  then  most  of  the  excited  atoms  will  deactivate 
before  reaching  the  lasing  cavity.  In  Chemical  Oxygen  Iodine  lasers  (COIL),  the  primary  flow  of  singlet  delta  oxygen 
(with  helium  carrier)  is  passed  through  a  supersonic  nozzle  to  the  laser  cavity.  The  fuel  (iodine  with  helium  carrier) 
is  injected  transversely  into  the  primary  flow  upstream  of  the  throat  in  subsonic  COIL  lasers  or  downstream  of  the 
throat  in  supersonic  COIL  lasers.  The  primary  advantage  of  supersonic  COIL  is  that  fuel  injection  is  decoupled  from 
the  upstream  singlet  oxygen  generator  (SOG).  Moreover,  fuel  and  carrier  injected  upstream  of  the  throat  in  subsonic 
COIL  must  pass  through  the  nozzle  throat.  Since  the  fuel  is  injected  downstream  of  the  throat  in  supersonic  COIL,  flow 
velocity  increases  and  the  partial  pressure  of  singlet  delta  oxygen  is  reduced,  both  of  which  reduce  the  rate  of  oxygen 
loss.  Several  designs  for  supersonic  COIL  lasers  have  been  described  in  the  open  literature  [5,  6,  7,  8,  9,  10,  11]. 
In  all  cases,  mixing  of  the  fuel  with  carrier  into  the  primary  flow  is  crucial  to  chemical  laser  design.  This  mixing  is 
achieved  by  injecting  one  of  the  reacting  species  into  the  main  flow.  In  turn,  this  problem  is  a  special  case  of  a  more 
general  problem  of  a  jet  in  crossflow.  The  latter  has  been  the  subject  of  some  detailed  experimental  studies  that  provide 
sufficient  data  for  experimental  validation  presented  in  the  following  sections  of  this  report. 

1.1  Shock  acceleration  of  a  gas  cylinder 

Richtmyer-Meshkov  instability  (RMI)  occurs  when  an  incident  shock  accelerates  an  interface  between  two  fluids  of 
different  densities  and  thus  amplifies  any  initial  perturbation  present  on  the  interface.  This  interfacial  instability  was 
theoretically  predicted  in  1960  by  a  Los  Alamos  theorist,  R.D.  Richtmyer  [12],  and  first  experimentally  observed  in 
1969  by  a  Russian  experimentalist,  E.E.  Meshkov  [13],  It  is  of  importance  to  applications  ranging  from  astrophysics 
[14],  inertial  confinement  fusion  [15]  to  supersonic  combustion  [16, 17],  The  fluid  configuration  leading  to  a  frequently 
considered  RMI  problem  is  shown  in  Fig.  1 .  Two  fluids  with  different  properties  are  separated  by  an  initially  diffuse 
interface.  The  shock  travels  from  light  fluid  (gray  color)  to  heavy  fluid  (black  color)  through  the  interface.  The 
development  of  RMI  for  this  configuration  can  be  divided  into  three  stages.  The  short  initial  stage  is  linear.  In  this 
stage,  the  incident  shock  wave  collides  with  a  perturbed  material  interface  and  bifurcates  into  a  transmitted  shock  and 
a  reflected  wave.  In  Fig.  1,  left,  the  pressure  and  density  gradients  are  locally  misaligned.  This  misalignment  leads 
to  baroclinic  generation  of  vorticity  at  the  interface.  The  resulting  vortex  roll-up  (Fig.  1,  right)  leads  to  growth  of  the 
perturbation  amplitude  of  the  interface.  The  perturbation  amplitude  is  very  small  as  compared  to  the  wavelength  of 
the  interface  and  hence  the  growth  of  interface  can  be  safely  assumed  to  be  linear  [12, 1 8],  Compressibility  effects  are 
important  in  this  stage  [12].  The  flow-field  in  this  stage  is  deterministic.  Later  in  this  stage,  spikes  and  bubbles  appear 
on  the  interface.  A  bubble  is  a  portion  of  the  light  fluid  penetrating  into  the  heavy  fluid  and  a  spike  is  a  portion  of  the 
heavy  fluid  penetrating  into  the  light  fluid. 
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The  second  stage  is  nonlinear  deterministic.  In  this  stage  the  spikes  and  bubbles  grow  substantially.  The  amplitude 
of  perturbation  grows  to  the  order  of  the  wavelength  and  hence  now  the  flow  is  nonlinear.  In  the  later  part  of  this  stage, 
roll-up  of  material  into  the  vortex  cores  occurs.  Roughly  at  the  same  time,  small  scales  also  appear  in  the  flow.  The 
small  scale  appearing  on  the  outer  edge  of  the  large  scale  vortex  core  may  be  due  to  shear-induced  secondary  instability 
while  the  small  scale  appearing  on  the  inner  edge  may  be  due  to  secondary  baroclinic  instability  [19,  20,  21].  In  this 
stage,  the  large  scales  in  the  flow  are  deterministic  [22,  23]  while  the  small  scales  are  stochastic.  Rightley  et  al.  [24] 
refer  to  this  stage  as  to  ordered  turbulence.  In  general,  CFD  codes  can  reliably  predict  the  large  scales  but  not  the 
small  scales  [25]. 

The  third  stage  is  disordered  turbulence  [25] .  In  this  stage,  the  secondary  instabilities  become  pronounced  and  lead 
to  the  onset  of  turbulence  with  chaotic  mixing.  Being  turbulent  in  nature,  this  stage  is  dominated  by  three-dimensional 
physics,  unlike  the  first  two  stages. 

The  development  of  RMI  depends  heavily  on  the  initial  conditions  that  influence  the  vorticity  deposition  in  the 
early  stages  of  RMI  evolution,  and  small  changes  in  the  initial  conditions  can  lead  to  profound  changes  in  the  late- time 
flow  morphology,  as  discussed  by  Budzinski  and  Benjamin  [26].  RMI  has  several  practical  applications  of  importance 
(some  of  which  were  mentioned  above),  including  strong  explosions  and  inertial  confinement  fusion,  which  makes 
reliable  numerical  simulation  of  RMI-driven  flows  a  priority.  At  the  same  time,  the  limited  energy  input  into  the 
RMI-driven  flow,  the  gradual  evolution  of  the  initially  deterministic  flow  towards  disorder  and  the  dependence  on  the 
initial  conditions  makes  RMI  an  attractive  problem  for  a  more  general  study  of  transition  to  turbulence.  Advances  in 
experimental  diagnostics  in  the  recent  years  have  produced  highly  resolved  quantitative  benchmarks  that  can  be  used  to 
validate  numerical  codes.  We  performed  a  detailed  comparison  of  a  2D  (two-dimensional)  numerical  simulation  of  an 
RMI-driven  flow  with  a  nominally  2D  initial  condition  geometry  commonly  used  for  code  validation  with  experimental 
results.  The  large-scale  flow  structure  observed  in  experiment  is  faithfully  reproduced  by  the  numerics.  The  numerical 
results  also  show  the  emergence  of  a  secondary  instability  which  has  been  observed  in  experiment.  However,  the 
statistics  of  the  small-scale  features  due  to  this  secondary  instability  are  quite  different  for  the  simulation  and  the 
experiment.  This  difference  is  likely  due  to  the  fact  that  the  scalings  for  the  disordered  features  are  dictated  by  the 
dimensionality  of  the  flow.  This  imposes  limitations  on  the  applicability  of  2D  simulations  to  flows  with  inherently 
3D  (three-dimensional)  small-scale  structure  evolving,  even  if  the  initial  conditions  and  the  large-scale  structure  of  the 
realizations  of  such  flows  are  effectively  2D. 

For  RMI  the  fluid  configuration  must  have  an  initially  perturbed  (non-planar)  density  interface  subjected  to  impul¬ 
sive  acceleration.  A  geometrically  simple  interface  (Fig.  1)  can  be  realized  experimentally  by  injecting  a  cylinder  of 
sulfur  hexafluoride  SF6  into  lighter  air.  The  interface  between  the  air  and  the  heavier  SF6  is  diffusive.  The  flow  is  im¬ 
pulsively  accelerated  by  a  planar  shock.  This  problem  has  been  the  subject  of  several  experimental  studies  [27, 28, 29], 
producing  flow  visualization  and  qualitative  data.  These  are  the  experimental  results  we  compare  with  the  results  of 
our  numerical  simulation. 

1.2  Supersonic  jet  in  crossflow  (JICF) 

Sonic  injection  into  a  crossflow  (whether  subsonic  or  supersonic)  that  occurs  in  chemical  lasers  is  a  special  case  of 
the  more  general  problem  of  a  Jet  in  Crossflow  (JICF)),  the  subject  of  the  current  work.  Mixing  is  expected  to  be 
dominated  by  vortical  structures  and  shock  structures  produced  by  this  interaction.  In  particular,  the  counter-rotating 
vortex  pair  (CRVP)  formed  in  the  jet-freestream  interaction  dominates  downstream  mixing.  Earlier  research  on  JICF 
was  focused  on  the  understanding  of  mixing  [30,  31,  32,  33,  34,  35].  Moreover,  reliable  experimental  data  for  JICF 
cases  is  available  for  comparison  with  numerical  simulations  of  a  flow  similar  to  that  which  occurs  in  a  chemical  laser 
nozzle.  Of  specific  interest  to  us  is  the  work  of  Gruber  et  al.  [35]  who  performed  an  experimental  mixing  study  of  a 
sonic  jet  injected  into  supersonic  crossflow  at  Mach  1 .98  revealing  the  overall  structure  of  the  flow.  The  main  flow  was 
comprised  of  air,  with  injection  either  of  air  or  of  helium  into  the  flow.  The  dominant  flow  structures  observed  in  a 
typical  JICF  configuration  involving  a  sonic  jet  and  supersonic  crossflow  are  shown  in  Figure  2.  The  adverse  pressure 
gradient  causes  the  turbulent  boundary  layer  to  separate  and  form  recirculating  zones  upstream  and  downstream  of 
injection.  The  jet  obstructs  the  supersonic  crossflow  and  thus  leads  to  the  formation  of  a  three-dimensional  bow  shock. 
The  underexpanded  sonic  jet  undergoes  expansion  forming  Prandtl-Meyer  expansion  fans  and  a  system  of  incident- 
reflected  oblique  shocks.  The  barrel  shock  is  the  barrel-shaped  structure  formed  due  to  the  expansion  fan  which 
terminates  in  the  Mach  disk.  The  horseshoe  vortex  and  wake  vortices  are  formed  downstream  of  the  injection.  The 
counter  rotating  pair  of  vortices  (CRVP)  dominates  the  cross-section  of  the  jet  in  the  far  field  and  improves  mixing  by 
entraining  the  crossflow  fluid.  The  dominant  structure  that  apprently  governs  the  mixing  in  both  validation  problems, 
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despite  many  differences  in  the  conditions,  is  a  counter-rotating  vortex  pair. 


2  Numerical  code 

GASP,  an  acronym  for  General  Aerodynamic  Simulation  Program,  is  a  CFD  code  from  AeroSoft  Inc.,  Blacksburg, 
VA.  GASP  is  a  structured  multi-block  CFD  solver  that  solves  an  integral  form  of  the  Reynolds  Averaged  Navier 
Stokes  (RANS)  equations  as  well  as  simplified  RANS  equations  such  as  thin-layer  RANS,  parabolized  RANS,  Euler 
equations  and  incompressible  RANS.  Density  and  pressure  are  time-averaged  (Reynolds-averaged)  while  temperature 
and  velocity  components  are  mass-weighted  averages  (or  Favre  averages).  For  compressible  or  reacting  flows  Favre 
averaging  is  favorable  since  it  incorporates  the  density  changes  while  averaging  the  Navier-Stokes  equations.  GASP 
uses  a  cell-centered  finite  volume  method.  GASP  has  a  variety  of  user-selectable  features  to  deal  with  compressible, 
mixing,  and  reacting  flows,  and  it  incorporates  several  turbulence  models  [37]. 

A  wide  variety  of  flows  has  been  studied  using  GASP  including  flows  in  a  combustion  chamber  and  flows  over 
reentry  vehicles.  GASP  code  includes  support  for  COIL  (chemical  oxygen  iodine  laser)  capability  developed  with 
support  from  AFRL  and  AFOSR.  GASP  is  being  used  by  NASA  for  verification  purposes  [38],  The  following  section 
describes  the  numerical  formulation  of  our  problems  with  GASP. 

3  Numerical  Formulation 

Our  simulations  involve  two  species  (air  and  SF6  for  the  first  problem,  air  and  helium  for  the  second),  and  the  flow  is 
considered  to  be  compressible  and  non-reacting.  As  mentioned  above,  we  use  a  commercial  CFD  solver  GASP  that 
solves  the  Navier-Stokes  equations  with  a  cell-centered  finite  volume  method.  The  fluxes  at  the  cell  faces  are  found 
using  the  MUSCL  (Monotonic  Upstream-centered  Scheme  for  Conservation  Laws)  scheme.  Several  flux  splitting 
methods  can  be  implemented  in  GASP,  with  the  possibility  to  assign  a  different  inviscid  flux  scheme  in  each  direction 
(x,  y,  z).  The  motivation  for  flux  splitting  is  provided  by  the  necessity  to  model  flows  with  local  subsonic  and  local 
supersonic  areas.  In  the  former,  flow  information  travels  both  upstream  and  downstream,  while  in  the  latter  it  is  down¬ 
stream  only,  which  can  pose  considerable  problems  unless  flux  splitting  is  implemented.  In  our  case,  the  interpolation 
of  primitive  variables  (from  cell  center  to  cell  faces)  is  third-order  upwind  biased  [37] .  In  the  case  of  the  first  validation 
problem,  the  flux  scheme  involves  an  assumption  of  no  flux  in  the  direction  normal  to  the  plane  of  the  flow  to  enforce 
two-dimensionality.  In  the  plane  of  the  flow,  Roe’s  upwind  inviscid  flux  scheme  [39]  is  used  for  the  first  validation 
problem.  Based  on  inviscid  2D  simulations,  the  Roe  with  Harten  flux  scheme  was  chosen  for  the  second  problem. 

The  viscous  terms  are  discretized  using  a  central  differencing  scheme  and  have  nominally  second-order  spatial 
accuracy.  The  laminar  transport  properties  (viscosity  and  thermal  conductivity)  are  specified  for  each  species.  For 
diffusivity,  Fick’s  law  of  binary  diffusion  is  used,  which  takes  into  account  only  the  effect  of  concentration  gradient. 
The  laminar  properties  are  calculated  for  each  species  in  the  flow.  Then  the  laminar  properties  for  the  mixture  are 
assessed  using  Wilkes’  mixing  model.  Species  are  assumed  to  have  different  densities,  but  pressure,  velocity  compo¬ 
nents  and  temperature  are  based  on  the  mixture  and  not  the  individual  species  {partial  multi-fluid formulation).  The 
binary  mixture  is  assumed  to  have  a  single  adiabatic  exponent. 

GASP  uses  either  an  explicit  Runge  Kutta  algorithm  or  an  implicit  dual  time  stepping  algorithm.  The  implicit 
dual  time  stepping  algorithm  has  two  steps.  In  the  first  step,  the  time  dependent  source  term  is  added  to  the  residual. 
During  the  second  step,  the  physical  time  is  kept  constant  while  the  residual  is  driven  to  zero.  GASP  uses  MPI  for 
parallelization.  The  computational  grid  is  imported  into  GASP  in  plot3d  format.  The  load  balancing  on  different 
processors  can  be  either  manual  or  automatic.  We  employed  the  dual  time-stepping  algorithm  with  automatic  load 
balancing. 

GASP  allows  the  use  of  several  turbulence  models,  from  algebraic  to  k  —  e,  k  —  u>  and  a  Reynolds  stress  model. 
In  the  case  of  a  shock-accelerated  flow  (validation  problem  1),  the  flow  is  initially  at  rest,  and  the  shock  interaction 
leads  to  formation  of  a  large-scale  vortex  structure  which  eventually  transitions  to  turbulence,  however,  small  scales 
take  time  to  develop.  Moreover,  from  experiment  we  have  a  priori  knowledge  [19,  27]  regarding  the  emergence  of 
small-scale  flow  features  in  our  transitional  flow,  and  in  the  time  interval  of  interest,  these  features  should  be  resolvable 
with  the  fine  mesh  described  later  in  this  report.  The  dimensionality  of  the  small  scales  is  also  discussed  later.  Thus  no 
turbulence  model  was  used  for  validation  problem  1 .  Viti  et  al.  [40]  found  that  for  a  problem  similar  to  our  second  JICF 
validation  problem  the  1998  k  -u>  turbulence  model  was  an  optimum  choice  considering  accuracy  and  computational 
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time.  Thus  we  selected  k  —  w  for  the  JICF  calculations.  The  wall  in  the  JICF  setup  was  assumed  smooth  by  setting 
the  dimensionless  roughness  to  Kj  —  1. 

The  Navier-Stokes  equations  describing  the  flow  under  these  assumptions  can  be  represented,  using  the  standard 
Einstein  notation,  as  shown  below. 


Mass  conservation  and  molecular  diffusion  : 


d  pn  ,  dpnUi  dpnVnt 

dt  dXi  dXi 

where  vni  is  mass  diffusion  velocity  of  the  species  n  in  i  direction.  Now  from  Fick’s  law  of  binary  diffusion  : 


PnVni  —  PmD 


nm 


dCn 

dxi 


where  Cn  represents  concentration  of  species  n  and  Dnm  is  the  diffusion  coefficient  of  species  n  with  respect  to  species 
m.  In  GASP  this  diffusion  coefficient  can  be  calculated  using  laminar  viscosity  (/q)  and  laminar  Schmidt  number  (Set) 
as : 


D 


nm  — 


nm 


We  used  a  constant  Set  =0.7  for  the  first  problem  (air  and  sulfur  hexafluoride).  For  the  second  problem,  the  laminar 
Schmidt  number  for  air  injection  was  0.7  (resulting  in  negligible  molecular  diffusion),  while  for  helium  injection  it 
was  0.3.  The  turbulent  Schmidt  number  was  0.5  for  all  cases.  The  simulation  was  performed  assuming  steady-state 
conditions.  The  pseudo-time  steps  for  the  second  problem  were  based  on  specifed  CFL  number  (local  at  first  and  then 
global). 

Momentum  conservation : 


dpui  d 
dt  +  dxi 


(pUiUi  +  pSij)  =  ~- 


(2) 


where  8tj  is  the  standard  Kronecker  delta  such  that : 


5ij  =  0  when  i  ^  j  and 


5ij  =  1  when  i  =  j 

and  Tij,  the  shear-stress,  given  by  (for  Newtonian  fluid) : 


Energy  Conservation  : 


where 


Tij  —  Hi 


f  dui  duj\ 
\dxj  ^  dxi ) 


2  duk 


dpe0  _5_ 
dt  dxi 


(pepUi  +pv,i  +  pe0Ui) 


d  ,  ,  dqt 

dxi  TyMj  dxi 


eo 

9i 


e  + 


UiUi 


(3) 


On  the  time  scale  of  interest,  the  contribution  of  heat  conduction  is  not  significant.  Thus  the  heat  transfer  coefficient  k 
was  set  at  a  constant  value  corresponding  to  that  of  air. 

Now,  since  we  have  two  species  (n,  m  :  1,2)  and  two  or  three  coordinates  (i,  j  :  1,  2  or  1,  2,  3),  the  equations 
(1),  (2)  and  (3)  can  be  finally  closed  by  assuming  each  of  the  species  to  be  thermally  and  calorically  perfect  and  using 
the  thermal  equation  of  state  : 

P  =  PRT  (4) 


along  with  the  thermodynamic  relation  : 


e  =  CvT 


(5) 
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By  using  the  Gauss  theorem,  the  above  system  of  equations  can  be  easily  represented  in  the  standard  integral 
form,  which  then  is  numerically  solved  in  GASP.  The  code  is  nominally  second-order  in  time  and  third-order  upwind- 
biased  in  space.  To  prevent  numerical  oscillations  near  the  shock  front,  we  used  a  Min-Mod  limiter  [37,  41].  The 
limiter  reduces  the  spatial  accuracy  to  first-order  in  the  vicinity  of  the  shock  and  hence  in  general  the  observed  spatial 
accuracy  may  be  lower  than  the  nominal  spatial  accuracy. 

4  Numerical  Simulation 

Details  specific  to  the  simulation  of  each  of  the  problems  we  considered  are  presented  below.  In  both  cases,  we  made 
an  effort  to  minimize  the  computational  resources  required  for  the  simulation. 

4.1  Shock  accelerated  gas  cylinder 

The  computational  domain  is  shown  in  Fig.  3.  The  domain  along  the  y-direction  was  reduced  to  1.0  cm  from  the 
physical  3.75  cm  half-width  of  the  shock  tube  by  applying  second-order  extrapolation  at  y  =  1.0  cm.  The  domain 
along  the  ^-direction  was  reduced  significantly  as  opposed  to  the  physical  length  of  the  shock  tube.  This  is  achieved 
by  using  the  upstream  condition  of  the  shock  in  the  test  section.  The  variation  of  pressure  in  the  shock  tube  is  shown 
in  Fig.  4. 

The  pressure  behind  the  shock  drops  from  p0  to  pi  and  hence  when  the  shock  hits  the  SFg  cylinder  in  the  test 
section,  it  has  upstream  conditions  corresponding  to  pressure  pi.  These  upstream  conditions  are  specified  as  the  initial 
conditions  for  the  shocked  region  (the  shaded  region  in  Fig.  3)  and  an  inflow  boundary  condition  is  used  at  x  =  0 
cm.  When  the  shock  hits  the  SF6  cylinder,  it  bifurcates  into  a  reflected  shock  wave  and  a  transmitted  shock  wave. 
The  reflected  shock  does  not  affect  the  development  of  RMI.  To  achieve  this  numerically,  we  change  the  boundary 
condition  at  x  =  0  to  reflective,  at  the  time  when  the  front  of  the  reflected  shock  wave  reaches  x  =  0  cm.  The  evolving 
SF6  cylinder  has  finite  velocity  and  to  ensure  that  it  stays  within  the  computational  domain,  a  uniform  negative  velocity 
is  superimposed  on  the  flow  in  the  x-direction.  It  is  important  to  note  that  the  SF6  cylinder  must  not  travel  towards  the 
x  =  0  boundary  because  that  would  lead  to  interaction  between  the  reflected  shock  and  SF6  cylinder. 

The  values  for  time  step  (At)  and  grid  spacing  (Ax)  were  found  in  a  convergence  study.  In  RMI,  the  moving  shock 
has  the  highest  propagating  speed  and  this  specifies  the  limiting  values  for  At  and  Ax.  Hence  a  one-dimensional 
moving  shock  problem  was  simulated  to  determine  At  and  Ax  .  Grid  sizes  of  100  pm,  200  pm  and  800  pm  were 
used.  After  obtaining  convergence  at  Ax  =  200  pm,  different  time  steps  (At)  were  used.  The  solution  converged  for 
At  =  0. 125  ps.  After  the  shock  passed  through  the  computational  domain,  the  value  for  At  was  increased  by  a  factor 
of  6. 


Table  1 :  Shock  tube  conditions 


Property 

Upstream 

Downstream 

of  shock 

of  shock 

Pressure  (Pa) 

121100 

80000 

x- velocity  (m/s) 

104.9 

0 

y- velocity  (m/s) 

0 

0 

Second-order  accuracy  in  time  and  third-order  accuracy  in  space  were  used  to  predict  RMI-driven  flow  evolution. 
As  it  was  mentioned  earlier,  the  effective  accuracy  of  the  code  is  lowered  by  the  use  of  limiters  in  the  immediate 
vicinity  of  the  shock.  However,  this  should  have  a  relatively  small  effect  on  the  evolution  of  the  flow  after  the  shock 
passage.  For  the  viscous  case,  the  Sutherland  viscosity  model  was  used  for  each  species,  and  the  mixture  kinematic 
viscosity  was  obtained  from  Wilkes’  law.  An  implicit  dual  time  stepping  Gauss-Seidel  scheme  was  employed  as  the 
temporal  algorithm.  The  initial  conditions  are  shown  in  Table  1 .  The  initial  distribution  of  SF6  is  based  on  a  simple 
Gaussian  fit  to  the  experimental  data  [42,  43],  with  the  peak  volume  fraction  of  SF6  at  80%. 
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Table  2:  JICF  conditions 

Quantity 

Case:  Air 

Case:  Helium 

Crossflow  Gas 

Air 

Air 

Injectant  Gas 

Air 

Helium 

Poo  (kPa) 

41.8 

41.8 

Too  (K) 

169 

168 

uQ o  (m/s) 

516 

515 

A/qq 

1.98 

1.98 

Poo  (kg/m3) 

0.860 

0.866 

Poo  (Pa-s) 

1.15  xlO-5 

1.14  xlO-5 

Reoo  -=  (PoodJoc 

>)/tio o(m_1)  3.87  xlO7 

3.91  xlO7 

Pj  (kPa) 

476 

405 

Tj  (K) 

250 

225 

Uj  (m/s) 

317 

882 

Mj 

1 

1 

Pj  (kg/m3) 

6.64 

0.867 

Pj  (Pa-s) 

1.60  xlO-5 

1.62  xl0~s 

Rej  =  (petj.Ue<j 

■Dj)/  pej  8.36  xlO5 

3.00  xlO5 

7 j 

1.4 

1.67 

j=^4- 

2.90 

2.90 

4.2  Jet  in  cross-flow 

The  conditions  for  our  simulations  were  selected  to  coincide  with  those  in  the  experiments  of  Gruber  et  al.  [35],  as 
summarized  in  Table  2  listing  the  static  values  of  pressure,  temperature,  etc.  for  the  two  cases  we  considered  (air  and 
helium  injection). 

The  computational  domain  is  shown  in  Fig.  5.  The  diameter  of  the  circular  injector  is  D  =  6.35  mm.  The  injector 
was  simulated  as  a  constant  area  duct  extending  from  —2  <=  <—  0.  The  structured  grid  was  generated  using 

Gridgen,  commercial  meshing  software  capable  of  generating  3D  structured  and  unstructured  grids  [44].  For  the  grid 
convergence  study,  we  used  D  =  6  mm  and  the  downstream  domain  extended  to  ^  =  10.  Our  fine  grid  had  4.7x10° 
grid  points.  Using  the  GASP  capability  of  removing  every  other  grid  point,  the  medium  grid  ( 6.1xlOr>  grid  points)  and 
coarse  grid  (S.lxlO4  grid  points)  were  generated  by  two  successive  coarsenings  of  the  fine  grid.  The  grid  had  about 
8  partially  independent  (in  one  direction)  zones,  with  the  total  number  of  zones  being  34.  Pressure  and  temperature 
were  monitored  at  the  wall  of  the  duct  (y  —  0,z  =  0)  to  determine  convergence  with  respect  to  the  number  of  cycles. 

A  no-slip  adiabatic  (split)  boundary  condition  was  imposed  at  the  wall  (-^  <=  0),  except  over  the  injection 
hole.  Turbulent  boundary  layer  profiles  of  primitive  variables  were  specified  as  crossflow  inflow  conditions  (at  fj  = 
—8).  These  boundary  layer  profiles  were  generated  using  the  space-marching  capability  of  GASP  and  were  verified 
with  boundary  layer  code  EDDYBL  [45],  To  match  the  experiments,  we  ensured  that  the  boundary  layer  thickness 
approaching  the  injector  was  about  ID.  Freestream  conditions  for  the  boundary  layer  were  the  same  as  those  listed  in 
Table  2.  Inlet  conditions  for  the  injector  flow  at  —  —  2  were  uniform  stagnation  pressure  and  temperature  across  the 
duct,  specified  to  produce  Mach  1  flow  at  the  injection  hole  =  0  and  to  match  the  experimental  static  pressure  and 
temperature.  The  no-slip  adiabatic  (split)  boundary  condition  was  applied  at  the  wall  of  the  injector  duct.  For  other 
outflow  boundary  conditions,  we  used  first-order  extrapolation. 

The  experiments  performed  by  Gruber  et  al.  [3  5]  and  Fric  and  Roshko  [46]  show  that  the  JICF  flow  field  is  highly 
unsteady  and  asymmetric.  In  particular  while  simulating  JICF,  Madden  and  Miller  [47]  found  that  the  assumption 
of  symmetry  and  steadiness  affects  the  development  of  wake  vortices  and  shear  layer  vortices.  However,  transient 
calculations  on  the  full  computational  domain  requires  prohibitive  computational  time  for  this  study.  Thus  we  assumed 
symmetry  and  steadiness  for  our  simulation. 
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5  Results 


The  results  from  simulations  are  qualitatively  and  quantitatively  compared  with  the  experiments  as  described  below 
for  each  problem. 

5.1  Shock-accelerated  gas  cylinder 

The  experimental  single-cylinder  data  from  the  LANL  DX-3  shock  tube  described  earlier  were  published  in  several 
papers  by  Vorobieff  et  a/. [27],  and  Tomkins  et  al.  [28,  29],  Here  we  mostly  use  the  work  of  Tomkins  et  al.  [29] 
for  validation.  Numerical  simulation  may  not  predict  the  physics  accurately  due  to  its  discrete  nature.  Also  the  2- 
dimensional  assumption  may  not  be  valid  at  late  times  when  secondary  instabilities  become  important.  Hence  we  have 
also  done  a  comparison  with  the  2-dimensional  simulations  done  by  Greenough  et  al.  [42]. 

The  widths  and  heights  of  the  mixing  zone  produced  by  shock  acceleration  of  the  SF6  cylinder  (Fig.  6a)  and  the 
spacing  between  the  vortex  cores  of  large  scale  counter-rotating  vortices  (Fig.  6b)  are  used  for  quantitative  comparison 
between  experiment  and  simulation.  Mesh  refinement  described  in  section  4.1  does  not  exert  much  influence  on  the 
dimensions  of  the  large-scale  vortices,  which  supports  the  use  of  the  one-dimensional  model  discussed  earlier  to 
determine  mesh  spacing.  However,  the  predictions  of  small-scale  secondary  instabilities  are  strongly  influenced  by 
mesh  spacing  (Fig.  7).  Note  that  the  secondary  instabilities  are  not  symmetric  and  hence  the  full  computational  domain 
must  be  used  unlike  the  one  shown  in  Fig.  3.  Decreasing  the  mesh  spacing  (Ax)  also  implies  decreasing  the  time 
step  (At)  so  as  to  keep  the  ratio  Ax/ At  constant.  Hence  when  the  mesh  spacing  is  halved,  the  computational  time 
increases  by  a  factor  of  8.  The  computations  were  performed  on  a  Linux  cluster  consisting  of  128  parallel  processors 
(256  MB  RAM  per  CPU).  The  computational  time  for  varying  mesh  spacings  is  shown  in  Table  3.  Considering  all 
these  factors,  we  decided  to  keep  the  final  mesh  spacing  as  50  /im.  Henceforth,  numerical  simulations  imply  a  full 
computational  domain  with  50  yam  as  the  mesh  spacing. 


Table  3:  Computational  Time  for  validation  problem  1. 


Mesh  Size  (yxm) 

Computational  Time 

200 

1.36  hours 

100 

10.9  hours 

50 

3.6  days 

25 

29  days 

12.5 

232  days 

The  qualitative  comparison  of  flow  fields  acquired  in  experiments  and  by  our  numerical  simulation  is  shown  in 
Fig.  8.  The  evolution  of  the  vortex  pair,  including  spikes  and  bubbles,  is  evident  in  our  simulations.  A  notch  at  the 
axis  of  symmetry  of  the  cylinder  appears  in  the  simulations  at  earlier  times,  as  can  be  seen  in  the  experiments  as  well. 
This  is  likely  due  to  shock  wave  focusing.  The  cylinder  appears  to  evolve  somewhat  slower  in  the  simulations  and  the 
small-scale  structures  appear  at  later  times  than  observed  in  experiments.  This  could  indicate  a  problem  with  vorticity 
deposition  during  the  initial  stage  of  the  RMI  simulation.  For  a  future  study  elucidating  the  issue,  it  would  be  highly 
desirable  to  procure  early-time  experimental  velocity  fields. 

The  secondary  shear  induced  instability  and  baroclinic  instability,  as  explained  earlier,  are  visible  in  our  results. 
The  secondary  instabilities  in  our  simulations  compare  well  with  these  from  the  simulations  of  Greenough  et  al.  [42] 
and  Rider  et  al.  [43]  using  the  Raptor  code  at  800  ya sec  with  comparable  mesh  sizes.  Raptor,  as  described  by  Howell 
and  Greenough  [48],  is  a  multi-physics  adaptive-mesh  refinement  (AMR)  code  that  uses  a  second-order  (both  in  space 
and  time)  Godunov  method  to  solve  Euler  equations. 

The  height  and  width  of  the  mixing  zone  (defined  as  the  horizontal  and  vertical  extent  of  the  perturbed  SF6  cylin¬ 
der),  and  the  vortex  spacing  are  used  for  quantitative  comparison  between  experiment,  our  simulation  and  the  sim¬ 
ulation  done  by  Greenough  et  al.  [42],  Vortex  spacing  from  experiments  was  not  available  for  early  times  and  was 
not  provided  by  Greenough  et  al.  [42],  These  comparisons  are  shown  in  Fig.  9.  The  results  show  that  the  overall 
trend  of  our  simulation  is  consistent  with  the  experiments  as  well  as  with  the  simulations  done  by  Grenough  et  al. 
[42],  There  is  some  difference  in  early- time  predictions  between  our  simulations  and  those  produced  by  the  Eulerian 
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AMR  codes,  with  the  latter  actually  reproducing  the  early-time  evolution  of  the  mixing  zone  width  better  (Fig.  9, 
top).  This  difference  could  be  due  to  the  limiter  used  in  our  code  to  avoid  fluctuations  near  the  shock  front.  At  later 
times,  however,  the  differences  between  our  results  and  the  experiments  are  very  close  to  those  between  predictions  of 
Greenough  et  al.  [42]  and  experiments.  This  illustrates  the  often-overlooked  point  that  quantitative  agreement  between 
two  simulations  does  not  necessarily  mean  that  the  same  degree  of  agreement  exists  between  either  of  the  simulations 
and  the  experiment. 

The  work  described  above  used  an  axisymmetric  Gaussian  density  distribution  for  the  initial  conditions.  However, 
in  experiments  the  initial  conditions  are  always  subjected  to  small  perturbations,  often  resulting  in  appreciable  devia¬ 
tions  from  radial  symmetry  (about  10  percent  in  terms  of  the  concentration  field,  according  to  Dwarkadas  et  al.  [49]). 
To  simulate  the  effect  of  such  perturbations  on  the  evolution  of  the  flow,  we  introduced  small  deviations  into  the  initial 
conditions  as  follows.  The  circular  distribution  was  replaced  with  a  slightly  elliptical  one,  with  the  ratio  between  the 
major  axis  and  the  minor  axis  1.05.  The  angle  9  between  the  major  axis  of  the  ellipse  and  the  direction  of  the  shock 
was  varied  in  increments  of  22.5°  between  0°  and  90°,  resulting  in  five  sets  of  initial  conditions. 

Figure  10  shows  the  evolution  of  the  flow  developing  from  these  initial  conditions  in  our  simulation.  At  early 
times,  variation  of  the  initial  conditions  has  little  influence  on  the  apparent  flow  structure.  The  only  flow  feature 
directly  affected  is  the  spike  produced  by  shock  focusing  (Fig.  10,  top  row  closeups).  Loss  of  symmetry  along  the 
centerline  in  the  direction  of  the  shock  propagation  leads  to  the  distortion  of  the  spike  very  similar  to  that  observed  in 
most  experiments. 

At  late  times,  the  overall  flow  morphology  remains  the  same,  manifesting  the  large-scale  structure  produced  by  the 
initial  RMI-induced  vorticity  deposition  and  the  smaller  structures  emerging  due  to  the  secondary  instabilities.  The 
specific  features  of  these  small-scale  structures  (e.g.,  position)  vary  for  different  6,  in  a  fashion  similar  to  that  observed 
in  experiment  for  initial  conditions  with  small  perturbations  inherently  present.  Thus  it  can  be  concluded  that  in 
the  “viewgraph  norm”  the  simulation  agrees  with  the  experiment.  Moreover,  the  simulation  captures  an  important 
feature  of  the  RMI-driven  cylinder  flow  in  the  stage  of  incipient  turbulence:  while  the  large-scale  structure  is  reliably 
reproducible  from  experiment  to  experiment,  the  small  scales  are  disordered.  In  the  simulation,  this  disordering  is 
produced  by  a  subtle  change  in  the  initial  conditions. 

In  the  analysis  of  experimental  measurements,  both  the  density  fields  [29]  and  the  velocity  fields  [27]  have  been 
ensemble-averaged.  One  could  expect  quantitative  agreement  between  these  ensemble-averaged  results  and  similar 
results  from  a  RANS  code.  By  subtracting  the  ensemble  average  from  the  instantaneous  experimental  fields,  fluctuation 
fields  were  produced  as  well.  Statistical  analysis  of  the  instantaneous,  ensemble-average  and  fluctuation  velocity  fields 
[27]  acquired  in  experiments  as  the  flow  transitions  to  turbulence  revealed  some  interesting  trends.  The  statistics  of 
the  flow  were  analyzed  in  terms  of  the  longitudinal  velocity  structure  functions 


Sn(r) 


(u(x  +  r)  -  u(x))  • 


where  u  is  the  velocity  vector,  r  =  |r|  and  <  . . .  >  denotes  ensemble-averaging.  For  the  second-order  velocity 
structure  function  S2  (r),  it  was  found  that  the  late- time  (~  1300  (is)  behavior  of  the  instantaneous  field  structure 
function  approaches  the  2/3  power  law,  consistent  with  Kolmogorov’s  1941  theory  of  fully  developed  turbulence. 
This  result  was  not  altogether  surprising,  as  the  flow  is  transitioning  to  turbulence.  What  was  more  interesting  was 
the  emergence  of  the  similar  behavior  at  a  much  earlier  time  (~  760  fis)  in  the  structure  functions  of  the  fluctuating 
velocity.  These  observations  were  based  on  the  measurements  of  two  velocity  components  in  a  plane  normal  to  the  axis 
of  the  cylinder,  with  an  effective  grid  resolution  of  about  70  fim.  The  question  that  arises  quite  naturally  is  whether 
the  same  trends  can  be  observed  in  the  analysis  of  our  numerical  data. 

Figure  1 1  shows  the  structure  functions  S2  (r)  produced  by  ensemble-averaging  of  the  velocity  fields  at  ~  900  /is 
developing  from  the  five  sets  of  elliptic  initial  conditions  described  above.  It  is  noteworthy  that  the  overall  agreement 
of  the  instantaneous  and  ensemble-averaged  structure  function  behavior  with  experiment  is  fairly  good.  Also  like  in 
the  experiment,  the  behavior  of  the  fluctuating  velocity  structure  function  is  distinctly  different.  Here,  however,  is 
where  the  similarity  between  experiment  and  theory  ends.  While  the  fluctuating  velocity  field  structure  function  in 
experiment  shows  a  trend  towards  2/3  power-law  scaling  starting  as  early  as  ~  760  /is,  the  fluctuating  components 
of  the  numerical  data  exhibit  a  distinctly  different  behavior,  with  the  fluctuation  structure  function  S2  growth  with 
r  falling  far  short  of  the  2/3  power  law.  Is  this  surprising?  Perhaps  it  should  not  be,  as  the  fluctuating  features  in 
experiments  are  three-dimensional,  while  similar  features  in  numerics  are  confined  to  two  dimensions,  and  no  simple 
increase  in  grid  resolution  is  likely  to  eliminate  this  problem,  even  if  the  2D  grid  resolves  the  Kolmogorov  scale.  In 
the  latter  case,  analysis  of  fluctuations  in  the  flow  developing  from  small  changes  in  the  initial  conditions  is  likely  to 
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reveal  statistics  consistent  with  those  predicted  and  experimentally  observed  for  spatially  two-dimensional  turbulence 
(see,  for  example,  Vorobieff  et  al.  [50]).  The  latter  leads  to  scalings  distinctly  different  from  those  predicted  by  the 
Kolmogorov  theory  and  observed  in  3D  experiments.  Thus  the  comparison  of  the  fluctuating  component  scalings  in 
experiments  and  numerics  reveals  an  important  physics  issue  of  limitations  inherently  present  in  a  2D  simulation  of  a 
spatially  3D  unstable  flow. 

In  simulations  performed  for  an  inviscid  case,  there  was  no  appreciable  difference  observed  in  the  results.  This 
supports  the  notion  that  viscous  dissipation  is  of  relatively  small  importance  on  the  time  scale  considered  here  (< 
1000  (is),  and  the  exact  choice  of  laminar  properties  does  not  have  much  influence  on  the  flow.  We  also  performed 
the  simulations  for  60%  peak  SF6  concentration.  In  this  case,  the  SF6  cylinder  seems  to  evolve  at  a  faster  rate  and 
the  secondary  instabilities  were  observed  at  times  much  earlier  than  the  80%  case.  This  can  be  attributed  to  steeper 
gradients  in  the  Gaussian  distribution  for  the  60%  case  (so  that  the  total  concentration  of  SF6  is  the  same  in  both  cases) 
as  compared  with  the  80%  case.  Steeper  gradients  in  density  fields  imply  stronger  baroclinic  deposition  and  hence 
faster  evolution  of  the  SF6  cylinder.  For  space  constraints,  the  results  are  not  presented  here.  The  peak  concentration 
of  SF6  in  the  initial  conditions  in  the  experiments  is  subject  to  some  fluctuations,  and  an  accurate  measurement  of 
the  initial  concentration  field,  as  well  as  the  early-time  velocity  field,  as  mentioned  earlier,  should  help  resolve  the 
discrepancy.  Zhang  et  al.  [21]  recently  adopted  a  feedback  approach  to  resolve  the  uncertainty  in  experiments. 
They  have  used  a  finite  thickness  layer  located  between  SF6  and  air.  The  parameters  of  this  finite  layer  are  changed 
with  iterations  until  good  agreement  is  achieved  between  the  simulations  and  experiments  [20,  21],  Zoldi  [36]  also 
concluded  such  need  to  change  the  initial  density  distribution  to  achieve  better  agreement  with  experiments.  However, 
such  approaches  cannot  be  used  to  validate  CFD  codes  as  they  rely  on  implicit  confidence  in  the  CFD  code. 

5.2  Jet  in  cross-flow 

For  the  JICF  problem,  we  also  performed  a  study  of  the  influence  of  the  mesh  size  on  the  characteristic  scales  of  the 
large-scale  flow  features  similar  to  the  one  described  in  the  beginning  of  section  5.1.  In  this  case,  the  characteristic 
scales  associated  with  the  counter-rotating  vortex  pair  dominating  the  flow  are  streamwise  and  spanwise  penetration 
(Fig.  12). 

Results  for  all  three  grids  described  in  section  4.2  show  similar  trends.  Transverse  penetration  is  nearly  grid 
independent  by  x/D  =  8,  while  spanwise  penetration  is  more  sensitive  to  grid  resolution.  At  x/D  =  9,  the  difference 
between  the  medium  and  fine  grid  results  is  just  3%.  The  area  enclosed  by  the  90%  injected  fluid  concentration  contour 
depends  on  the  shape  of  the  jet  in  the  CRVP  as  well  as  the  penetration  quantities  that  indicate  its  outer  dimensions. 
There  is  signifcant  numerical  diffusion  on  the  coarse  grid  that  is  manifested  in  larger  spanwise  penetration  and  enclosed 
area.  Thus  it  is  clear  that  the  coarse  grid  does  not  resolve  the  jet  cross-section  adequately.  The  medium  grid  was  found 
to  resolve  the  shock  and  vortical  structures  satisfactorily  and  was  used  for  the  remainder  of  the  simulations.  Grid 
convergence  computations  show  that  the  order  of  convergence  is  just  under  2,  in  agreement  with  more  extensive  grid 
convergence  studies  using  GASR 

The  flow  features  shown  in  Fig.  2  are  adequately  resolved  in  our  simulations.  These  flow  features  are  quite  similar 
for  the  two  injected  fluids  (air  and  helium).  Figure  13  shows  Mach  number  contours  for  air  and  helium  injection,  with 
which  the  upstream  separation  region,  barrel  shock,  bow  shock  and  Mach  disk  identified.  Helium-injected  flow  has  a 
stronger  expansion  in  the  barrel  shock,  reaching  about  Mach  7,  while  the  flow  with  air  injection  reaches  about  Mach 
6.  For  helium,  the  Mach  disk  is  centered  at  x/D  «  1.8,  while  for  air  it  is  centered  at  x/D  «  2.2.  Apart  from  these 
differences,  the  flow  structure  of  the  two  cases  is  very  similar. 

In  both  cases,  a  horseshoe  vortex  system  forms  upstream  of  the  jet,  as  shown  by  streamtrace  visualization  in  Fig.  14. 
An  (imaginary)  cylinder  of  diameter  D  and  length  1.2D  located  at  the  origin  is  shown  for  reference.  The  horseshoe 
vortex  system  originates  in  the  lower  part  of  the  upstream  separated  boundary  layer.  The  crossflow  is  effectively 
blocked  by  the  presence  of  the  jet  so  that  incoming  fluid  enters  a  recirculating  region  and  then  passes  to  one  side  of 
the  jet.  Upstream  of  the  jet  near  the  symmetry  plane,  the  horseshoe  vortex  has  negative  z  vorticity  which  is  turned  in 
the  streamwise  direction  by  the  crossflow.  While  Fig.  14  shows  the  streamtraces  for  air  injection,  the  horseshoe  vortex 
for  helium  injection  is  veiy  similar. 

The  wake  vortex  system[46]  for  the  air  injection  case  is  shown  by  streamtraces  in  Fig.  1 5  along  with  the  (imaginary) 
reference  cylinder.  The  origin  of  the  wake  vortex  is  in  another  recirculating  zone  in  the  upper  part  of  the  upstream 
separated  boundary  layer.  Upstream  of  the  jet  near  the  plane  of  symmetry,  the  wake  vortex  has  positive  a  vorticity. 
The  wake  vortex  wraps  tightly  around  the  jet  to  converge  in  its  wake  at  x/D  «  3,  and  then  each  leg  diverges  in  the 
spanwise  direction  further  downstream.  Thus  the  spanwise  vorticity  is  quickly  converted  to  streamwise  vorticity.  It 
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should  be  noted  that  there  is  a  single  pair  of  wake  vortices  downstream  of  the  jet.  This  is  not  obvious  from  the  figure 
based  on  selected  streamtraces.  Although  the  wake  vortex  originates  well  above  the  horseshoe  vortex,  its  legs  approach 
the  wall  downstream  of  the  jet,  while  the  horseshoe  vortex  system  remains  at  a  nearly  constant  vertical  height. 

Results  of  the  experiments  and  our  predictions  may  be  compared  so  long  as  it  is  noted  that  the  experimental  images 
indicate  intensity  of  the  Rayleigh/Mie  scattering  from  particles  in  the  crossflow  (with  black  to  white  reversal).  Gruber 
et  al.  [30,  35]  hypothesize  that  the  intensity  is  directly  proportional  to  the  concentration  of  particles  in  the  mixed 
fluid,  and  is  thus  proportional  to  the  mass  fraction  of  crossflow  fluid.  With  the  black  to  white  reversal,  white  in  the 
images  corresponds  to  mass  fraction  of  jet  fluid.  No  quantitative  relationship  between  intensity  of  scattered  light  and 
particle  concentration  was  reported.  With  this  taken  into  consideration,  our  time-averaged  predictions  of  crossflow 
mass  fraction  are  compared  with  ensemble-averaged  experimental  images. 

Figure  16  compares  images  from  the  experiments  and  our  predictions  of  mass  fraction  of  the  crossflow  fluid  at 
several  locations  in  the  streamwise  direction  for  both  air  and  helium  injections.  Unlike  the  simulations,  the  experi¬ 
mental  images  do  not  resolve  the  fine  structure  associated  with  the  jet  cross-section  in  these  end  views  of  the  y-z  plane. 
Our  results  show  good  agreement  with  the  experiments  for  the  large-scale  structure  of  the  jet.  The  counter-rotating 
vortex  pair  (CRVP)  becomes  evident  by  x/D  =  4.  From  our  results,  the  case  of  helium  injection  has  smaller  regions 
of  unmixed  jet  fluid  as  compared  to  air  injection.  At  x/D  =  4,  helium  injection  also  shows  the  significant  effect  of  a 
secondary  instability  (with  smaller  scale)  upon  the  primary  CRVP  instability.  This  is  evident  to  a  lesser  extent  for  air. 
This  secondary  instability  grows  more  rapidly  downstream  for  helium  injection,  although  the  shape  of  the  jet  cross- 
section  is  similar  between  the  two  cases.  Quantitative  comparison  between  the  penetration  parameters  in  experiment 
and  numerics  (Fig.  17)  shows  a  fairly  good  agreement,  especially  taking  into  consideration  the  rather  limited  nature  of 
the  experimental  data. 

6  Conclusion 

We  have  performed  two-dimensional  numerical  simulations  of  two  validation  problems  involving  shock  acceleration 
of  an  initially  diffuse  heavy-gas  (SF6)  cylinder  embedded  in  lighter  gas  (air),  and  injection  of  a  supersonic  jet  of 
gas  (air  or  helium)  into  supersonic  crossflow  (air).  Unlike  many  earlier  simulations,  this  work  uses  a  robust  and 
hydrodynamically  relatively  simple  code  (e.g. ,  no  adaptive  mesh  refinement)  and  limited  computational  resources, 
nevertheless  producing  overall  good  agreement  with  experiment. 

The  main  conclusion  of  the  validation  exercises  is  that  the  hydrodynamic  model  of  GASP  is  adequate  for  the 
numerical  simulation  of  high-speed  compressible  mixing  flows  in  chemical  lasers.  Two-dimensional  simulations  of 
intrinsically  three-dimensional  flow  phenomena  might  be  useful  in  predicting  the  macroscopic  flow  features,  however, 
they  cannot  be  expected  to  give  a  correct  statistical  prediction  of  the  small  scales  important  for  mixing.  In  general, 
it  is  reasonalbe  to  expect  the  geometric  parameters  of  large-scale  flow  features  (counter-rotating  vortex  pairs)  to  be 
predicted  with  an  error  not  exceeding  10-15%.  The  grid  parameters  and  computational  setup  described  here  for  the  jet 
in  crossflow  should  be  adaptable  to  a  specific  chemical  laser  injector  geometry. 
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Figure  1 :  Schematic  representation  of  a  planar  shock  acceleration  of  an  initially  diffuse  cylindrical  density  interface: 
left  -  prior  to  acceleration,  right  -  after  acceleration.  Note  the  initial  misalignment  between  gradients  of  pressure 
(black  arrows)  and  density  (gray  arrows)  leading  to  deposition  of  vorticity  (white  dashed  arrows).  Note  that  this 
schematic  does  not  represent  many  of  the  subtle  effects  in  the  flow,  e.g.,  shock  reflection  or  focusing. 


Figure  2:  Schematic  diagrams  showing  the  vortical  and  shock  structure  of  the  JICF  flow  field  (from  Ref.  [35]):  (a) 
side  view,  (b)  perspective  view. 
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Figure  3:  Computational  domain  for  validation  problem  1  (shock-accelerated  gas  cylinder).  Note  that  only  the  top  half 
of  the  actual  domain  is  shown. 


Figure  4:  Pressure  variation  along  the  shock  tube. 
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Figure  5:  Schematic  of  the  computational  domain  for  validation  problem  2  (jet  in  cross-flow)  showing  various  grid 
zones. 
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Figure  6:  (a)  The  width  and  height  of  the  mixing  zone  produced  by  shock-acceleration  of  the  SF6  cylinder,  (b)  Vortex 
spacing.  Dark  region  corresponds  to  air  in  density  field  and  low  pressure  in  pressure  field. 
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Figure  8:  Comparison  of  experiment  (top)  and  our  numerical  simulation  (bottom).  In  experimental  images  dark  color 
corresponds  to  higher  SF6  concentration.  In  numerics  blue  color  corresponds  to  air. 
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Figure  9:  Comparison  of  mixing  zone  height  (top),  width  (middle)  and  vortex  spacing  (bottom)  between  experiment 
and  simulation. 
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Figure  10:  Flow  morphology  for  different  initial  conditions.  Top  :  at  time  =  345  (is  and  bottom :  at  time  =  900  (is.  The 
first  image  in  both  rows  corresponds  to  the  circular  distributions  while  all  other  images  are  for  elliptical  conditions. 
Closeup  shows  the  notch. 


Figure  1 1 :  Second-order  longitudinal  velocity  structure  functions  of  instantaneous,  ensemble-averaged  and  fluctuating 
velocity  fields  obtained  from  elliptic  ICs  (see  text)  at  ~  900  (is.  The  structure  functions  are  normalized  by  their 
average  value,  the  scale  r  is  normalized  by  the  nominal  gas  cylinder  diameter  Dq.  The  second-order  velocity  structure 
function  obtained  in  experiment  [27]  at  ~  760  (is  is  shown  for  comparison.  Experimental  and  numerical  data  are 
labeled  ‘E’  and  ‘N’  respectively  in  the  captions. 
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Figure  12:  Schematic  of  an  arbitrary  end  view  (at  fixed  downstream  location  x/D )  showing  crossflow  fluid  mass 
fraction  contours  and  transverse  and  spanwise  penetration  scales.  Red  denotes  pure  jet  fluid  and  blue  denotes  pure 
crossflow  fluid.  D  is  the  diameter  of  the  injector. 
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Figure  13:  Jet  in  crossflow:  Mach  number  contours  in  the  x  —  y  plane  at  the  centerline  (z  —  0,  or  plane  of  symmetry) 
for  two  injected  fluids  (air,  left,  and  helium,  right).  Shock  structures  and  upstream  recirculating  region  are  identified. 
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Figure  14:  Jet  in  crossflow  (air  injection):  streamtraces  showing  the  horseshoe  vortex  system,  along  with  an  (imagi¬ 
nary)  cylinder  of  diameter  D  and  length  1 .2  D  located  at  the  origin  shown  for  reference. 
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Figure  16:  Mass  fraction  of  the  crossflow  fluid  for  end  views  ( y  —  z  plane)  at  downstream  distances  x/D  =  4  (left)  and 
x/D  =  8  (right).  The  left  column  shows  experimental  measurements  of  Gruber  et  al.  [30,  35]  with  our  predictions  in 
the  right  column.  Black  corresponds  to  100%  crossflow  fluid,  while  white  corresponds  to  0%  crossflow  fluid  (or  100% 
jet  fluid).  For  each  downstream  position,  results  for  air  injection  are  shown  on  top  and  results  for  helium  injection  are 
shown  on  the  bottom. 


Figure  17:  Jet  in  crossflow  penetration  based  on  a  side  view  of  the  x  —  y  plane  at  z  =  0.  Here  the  jet-to-crossflow 
momentum  ratio  is  J  =  2.90  and  a  =  D/2. 
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